Quantum phases of hardcore bosons with long-range interactions on a square lattice 
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We study the ground-state phase diagrams of hardcore bosons with long-range interactions on 
a square lattice using the linear spin-wave theory and a cluster mean-field method. Specifically, 
we consider the two types of long-range interaction: One consists only of the nearest- and next- 
nearest-neighbor interactions, and the other is the dipole-dipole interaction that decays with the 
interparticle distance r as ~ r~ 3 . It is known from previous analyses by quantum Monte Carlo 
methods that a checkerboard supersolid (CSS) is absent in the ground-state phase diagram of the 
former case while it is present in the latter. In the former, we find that quantum fluctuations around 
mean-field solutions are enhanced by the direct competition between the checkerboard and striped 
solid orders and that they destabilize the CSS phase. On the other hand, the emergence of the 
CSS phase in the latter case can be attributed to the absence of such a competition with other 
solid orders. We also show that the cluster mean-field method allows for the determination of phase 
boundaries in a precise quantitative manner when scaling with respect to the cluster size is taken 
into account. It is found that the phase transition between the superfluid and the solid (or CSS) is 
of the first order in the vicinity of the particle-hole symmetric line. 

PACS numbers: 03.75.-b, 05.30.Jp, 67.80.kb 
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I. INTRODUCTION 



Can a solid exhibit superfluidity in lattice systems? 
This question was first investigated theoretically by Mat- 
suda and Tsuneto^ in the context of the quantum 
lattice-gas model for 4 He, which assumes that atoms 
move only on fixed lattice points even in the liquid 
phased Using the lattice representation, they discussed 
the possibility of supcrsolidity, which is characterized by 
the coexistence of solid (diagonal) and superfluid (off- 
diagonal) long-range orders, in bulk and thin film of 4 Hc. 
In the lattice system, the continuous translational in- 
variance of the system is broken by the presence of the 
background discrete structure, and the "solid" means a 
state in which a discrete translational invariance is bro- 
ken spontaneously. Recently, this issue has attracted re- 
newed interest in connection with ultracold Bose gases 
in optical lattices. The creation of gases with strong 
dipole-dipole interactions 4 ^— has provided an ingredient 
essential for the emergence of supersolid phases, namely 
long-range interactions. Moreover, the precise control- 
lability of optical-lattice systems has inspired theoreti- 
cal explorations of supersolid phases in various types of 
lattice structure, such as chain ) 10 ! 11 square^ - — trian- 
gular; 2 ^— honeycomb ) 30 ' 31 kagome^ and cubi o 16 ' 33 " — 
lattices. We also note that the formation of checkerboard 
density-wave order has been experimentally observed in 
Bose-Einstein condensates coupled with an optical cav- 
ity^ 

For understanding lattice supersolids, it is important 
to address the following questions: in what situations 
the coexistent state can emerge and why it can be stable 
in such situations. Extensive studies over the past few 
decades have provided answers to these questions. For 



example, previous researches demonstrated that no su- 
persolid phases can exist in the ground-state phase di- 
agram of the hardcore Bose-Hubbard model with the 
nearest-neighbor (NN) interaction for bipartite lattices 
such as square 3 - 7 - and honeycomb 30 ' 31 lattices. In these 
cases, uniform supersolid states are unstable towards the 
formation of domain walls ) 14 ' 38 and the system under- 
goes phase separation into superfluid and solid phases. 
In order for supersolid phases to be present, one has to 
modify the model by, e.g., introducing dipole-dipole in- 
teractions^ 3 ^ or treating softcore bosons . 13 ' 14 In con- 
trast, the triangular-lattice system of hardcore bosons 
with only the NN interaction has stable supersolid phases 
for the fillings 1/3 < p < 2/2,M^r2^ As for the case of 
the kagome lattice, although the mean-field (MF) anal- 
ysis predicts the existence of supersolid states^ - they 
are destabilized by the effects of strong quantum fluc- 
tuations.— 

In this paper, focusing on the supersolid phase with 
checkerboard solid order, we analyze ground-state prop- 
erties of hardcore bosons with long-range interactions on 
a square lattice by means of the linear spin-wave (LSW) 
theory and a cluster mean-field (CMF) method. In this 
system, the range of the interactions makes a qualita- 
tive difference in the emergence of checkerboard super- 
solid (CSS) states. The previous quantum Monte Carlo 
(QMC) calculations 3 - 7 - have shown that no CSS phase is 
present between the superfluid (SF) and checkerboard 
solid (CS) phases in the system with only the NN inter- 
action Ui and the next-nearest-neighbor (NNN) interac- 
tion V2 ■ On the other hand, it is known that the infinite- 
range dipole-dipole interaction, which decays as the in- 
verse cube of the distance, can stabilize the CSS states.— 
We will clarify the reasons why the dipole-dipole interac- 
tion can stabilize the CSS states unlike the case of only 
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the NN and NNN interactions. 

The MF ground-state (classical) properties of the hard- 
core Bose-Hubbard models with dipole-dipole interaction 
and with only the NN and NNN interactions have already 
been discussed separately in previous works] 19 ' 41 " — We 
first review those results from the standpoint of com- 
paring the two types of interactions. When assuming 
that the system is in the phases with checkerboard (two- 
sublattice) order, the MF energy of the dipolar model 
can be naturally written in the same form as that of 
the model with effective NN and NNN interactions, V* 
and V£ s . We find that the value of V£ s /V£ s is very 
large, and it leads to a large region of CSS phase in the 
ground-state phase diagram at the MF level. Second, 
from the LSW analysis, we show that quantum fluctu- 
ations around the MF solutions are not so strong com- 
pared to the case of only the NN and NNN interactions, 
which is attributed to the absence of the direct compe- 
tition between the checkerboard and other solid orders. 
These two factors lead to the emergence of the stable 
CSS state in the dipolar system unlike the case of the 
shorter-range interactions. 

Moreover, including the effects of quantum fluctua- 
tions, we derive the ground-state phase diagrams. Al- 
though some of the results have already been known from 
previous QMC works, we reconsider the issue in detail in 
terms of another numerical approach based on a large- 
size CMF method^ From a comparison with the QMC 
result 1 ^ for the model with dipole-dipole interaction, it is 
shown that the CMF method combined with cluster-size 
scaling can locate the phase boundaries quantitatively. 
We also derive the phase diagram of the model with only 
the NN and NNN interactions and confirm that the re- 
gion of stable CSS phase almost completely disappears 
due to the strong quantum fluctuations. Moreover, we 
find the first-order phase transition between the SF and 
the CS (or the CSS) in the close vicinity of the particle- 
hole symmetry line for the both models. It is worth 
stressing that our CMF procedure is free from the minus- 
sign problem even when applying to frustrated systems. 
Moreover, it is useful to study metastability phenom- 
ena such as hysteresis^ since one can get all station- 
ary points of the free energy including metastable and 
saddle-point solutions. 

The remainder of the paper is organized as follows. 
In Sec. UH we introduce our models describing hardcore 
bosons with two types of long-range interactions in a 
square lattice. In Sec. Hill we show the ground-state 
phase diagrams of the two models within the mean-field 
theory. In Sec. IIV1 we perform the LSW analyses to dis- 
cuss the strength of quantum fluctuations around the MF 
solutions. In Sec. [V] applying a CMF method and the 
cluster-size scaling, we obtain the phase diagrams includ- 
ing the effects of the quantum fluctuations. Moreover, we 
summarize the reasons why the dipole-dipole interaction 
stabilizes the CSS states, based on the results obtained 
in Sees. IHUfVl The conclusion is given in Sec. IVI1 
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FIG. 1: Schematic pictures of the sublattice structure of the 
(I) checkerboard and (II) stripe patterns. The circles indicate 
the sites of the square lattice and the lattice sites of the same 
color belong to the same sublattice. 



II. HARDCORE BOSE-HUBBARD MODELS 

We consider interacting hardcore bosons on a square 
lattice given by the following Hamiltonian: 



H 



-J Yl^)® 1 



(1) 



U,i) 



3,1 



where a- and 



operators of the hardcore bosons at site j, J denotes the 
hopping amplitude between NN pairs, and fi the chemical 
potential. The hardcore boson limit means the situation 
where two or more bosons are not allowed to occupy the 
same site due to the strong on-site interaction U — > oo. 
We assume the existence of a long-range interaction Vij 
between the hardcore bosons and consider two different 
forms of Vji such that we study the effect of long-range 
interactions on the stability of supersolid phases through 
the comparison of the two models. 
The first one is given by 

r vx {\v j -v l \ = d) 

Vjl = < V 2 (\v J -v l \ = ^2d) (2) 
I (otherwise) [Vi — V<i model]. 

where d is the lattice spacing and Tj = (j x d,j y d) with 
integers j x and j y is a lattice vector at site j. The pa- 
rameters Vi > and V2 > represent the strength of the 
NN and NNN interactions, respectively. The NN inter- 
action V\ tends to induce the checkerboard density-wave 
order depicted in Fig. [Tjl), while the strong NNN in- 
teraction V2 favors the stripe pattern in Fig. [TJII) i 37 i 43 
Thus, the Hamiltonian in Eq. (fTJ) with Eq. ((2|), which we 
refer to as the "V1-V2 model," is a minimal model for 
studying the competition among two different solid or- 
ders and superfluidity induced by the hopping J . 37 i 41 ~ — 
We will focus on the regime of checkerboard ordering, 
V-xjVx < 1/2^-ii. 

The second one is the isotropic dipole-dipole interac- 
tion that is more realistic from an experimental point 
of view. In experiments of ultracold gases, one of the 
most promising way to prepare long-range interacting 
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systems is the use of the so-cold "dipolar" atoms, such 
as chromium^ dysprosium^ and erbium^ or molecules, 
such as KRb£>& and LiCs.— These atoms and molecules 
have a large (magnetic or electric) dipole moment, which 
leads to strong long-range forces among the dipolar par- 
ticles. We assume that the dipole moments are fully po- 
larized along the direction perpendicular to the lattice 
plane. In this case, the interaction between the dipolcs 
works isotropically and its long-range part can be well 
approximated by 



Vji 



Vd 3 / \rj 




(.1 * 

U = 



[V dip model]. (3) 



We refer to the model given by the Hamiltonian in Eq. JT]) 
with Eq. ([3]) as the "Vdi p model," hereafter. The dipolc- 
dipole interaction falls off as the inverse cube of the 
distance as {V, 0.354V, 0.125V, 0.089V, • • •} (for the NN, 
NNN, third, fourth neighbors). Therefore, it appears 
that most of the essential physics can be captured with 
just the first two terms, namely within the V1-V2 model. 
In fact, as will be shown in the next section, the MF 
phase diagrams of the two models are very similar; there 
are regions of the standard SF phase, solid phases, and 
supersolid phases, including the CSS phase. 

However, previous numerical analyses based on the 
QMC method demonstrated that the correct ground- 
state phase diagrams, which include quantum fluctua- 
tions, have a crucial difference between the finite-range 
V1-V2 model and the infinite-ran ge V dip model. For the 
V1-V2 model, the authors of Ref. |37| concluded that the 
CSS phase predicted by the MF theory is completely 
destabilized by strong quantum fluctuations and does not 
appear in the QMC calculations for any value of V2/V1, 
although they checked it only for Vi = 3 J. In contrast, as 
for the Vdip model, the main features of the MF phase di- 
agram can survive^ including the existence of the CSS 
phase. This indicates that the long-range part of the 
dipole-dipole interactions plays a crucial role in stabiliz- 
ing the CSS phase. In the following sections, we shall 
analyze the ground states of the two models and discuss 
the role of the long-range interactions in the emergence 
of the CSS state in order to clarify the reasons why the 
two models have the qualitative difference. 

The instability of supersolid phases against phase sepa- 
ration has often been discussed from a perturbative point 
of view assuming that J/Vji is small; the total energy 
gains from the lowest-order hopping process of doped 
bosons (or holes) and from the surface energy are com- 
pared on classical solid states with/without a domain- 
W all,Mi22i2Si22 However, in the V dip model, the CSS states 
can appear even for relatively large values of J/Vji, and 
the structure called the devil's staircase with many dif- 
ferent types of solid states emerges in the region of small 
values of J/VjiJ£ Hence, we will present a more careful 
discussion from a different angle by using the LSW and 
CMF methods. 



III. CLASSICAL GROUND STATES 

To begin with, we show the ground-state properties 
within the MF theory. From the equivalence of the 
hardcore-boson and spin-1/2 operators^ Eq. (fT]) can be 
mapped onto the spin-1/2 XX Z model with long-range 
Ising-type interactions: 



-2J 



2 ^ 

3,1 



(4) 



j — (SJ, Sj, Sj) is the pseudospin operator which 



where S 

satisfies the commutation relations 



[S?,Sf] =ic M „A,S#fy. (5) 

In the pseudospin language, the occupied and unoc- 
cupied states of bosons correspond to the spin-up and 
spin-down states, respectively. Thus the filling factor, 
which is the average density per site, of hardcore bosons 
can be calculated through the relation p = ^2j(fij) /M = 

1/2 + Y,j(Sj)/ M - Hcrc > M is thc number of lattice 
sites. The pseudospin raising and lowering operators 
Sj+iSj play the role of the creation and an- 



nihilation of the hardcore bosons; a] 



SI 



Gj = S J 



3 3 

The effective magnetic field acting on the pseudospins is 



given by h = (J, — zV /2 with 



Vi + V 2 (for the V1-V2 model) 
2.258V (for the V dip model) 



(6) 



where z = 4 is the coordination number of the square 
lattice. As can be obviously seen from the definition, the 
zero magnetic field corresponds to the particle-hole sym- 
metric point (/i = zV/2) of the hardcore-boson model. 
Moreover, the density pj and the condensate wave func- 
tion &j of bosons are expressed by the longitudinal com- 
ponents (Sj) and the transverse components (S^~) as 

pj = +1/2 and *j = From these corre- 

spondences, we can use the calculation methods which 
have been developed in the field of quantum spins for 
studying hardcore-boson systems. 

At zero temperature, replacing the local pseudospin 
operators in Eq. (|4]) with the classical vectors of length 
S = 1/2, 



QCl 



S^cos ipj sin 9j , sin ifj sin 0j , cos 6j), (7) 



we obtain the MF (classical) energy as a function of the 
orientation of the local pseudospins {8j, 

S 2 

3,1 



s 2 r 

E = -— ^2 \2Jji smOj sin0; cos(<pj - </?/) 



-Vji cos 6j cos 9 1 



hS ^ cos 9j , 



(8) 



4 



where 

J | (otherwise) v ' 

Minimizing the MF energy with respect to {9j, (fj}, we 
derive the classical pseudospin configurations in the usual 
manner ^ir— and translate the results into the hardcore- 
boson language. Without loss of generality, we can take 
(fj = 0, which means that the canted spins are assumed 
to lie in the xz plane. The procedure described here 
gives the same results as those obtained by the standard 
decoupling technique for the intersite spin-exchange in- 
teraction terms (i.e., the Weiss molecular- field theory) at 
T = 0. 



A. The MF results for the V1-V2 model 

In this subsection, let us briefly review the MF results 
for the V1-V2 model4i~— We mainly focus on the phases 
with the two-sublattice structure described in Fig. [HI). 
Within this checkerboard structure, we can describe the 
CS and CSS states in addition to the uniform SF state. 
The CS state, which is an insulating state appearing 
at half filling, has the checkerboard density- wave or- 
der characterized by pq = Y^ji^j) °xp(iQ • rj)/M with 
Q = (Tr/d,Tr/d), while the SF state has the off-diagonal 
long-range order characterized by = /M ■ The 

CSS state has both of the two (diagonal and off-diagonal) 
orders. In addition, completely empty (p = 0) and fully 
occupied (p = 1) states also appear. These trivial in- 
compressible states can be regarded as a kind of Mott 
insulator (MI) states. 

In the classical limit, these states can be expressed in 
terms of pseudospin angles as follows: 

cos6» A = -cos6» B = 1 (CS) 

6 A ^ 9 B and sin6»A,sin6»B + (CSS) , 

sm6 A =sin0 B ^0 (SF) ' [ } 

cos 9 A = cos 6 B = 1 or - 1 (MI) 

where 6 a and 8b are the canting angles of the pseudospins 
on sublattices A and B [see Fig. Hjl)] . The MF energy in 
Eq. ((5J) per site can be rewritten as a function of 9a and 
9b- 

E^ ch) /M = -4JS 2 sin 9 a sin 9 B + 2Vi S 2 cos 9 A cos 9 B 
+V 2 S 2 (cos 2 9 A +cos 2 9 B ) 
-hS(cos9 A + cos9 B )/2, (11) 

and the filling factor is given by p = 1/2 + S (cos 9 a + 
cos$b)/2. The ground-state phases are determined so as 
to minimize the MF energy Eq' with respect to 9 a and 
9b- Figures [2j a) and HJb) show the phase diagrams in 
the (J/Vi, /i/Vi)-plane for two different values of V2/V1. 
The value of J at the tip of the CS phase, J c , is given by 
(Vi — V^)/2 within the MF theory. The phase boundaries 



between CS and CSS, between CSS and SF, and between 
SF and MI are given by ft. = ±/i c i, ih C 2, and ±/i C 3, where 

h cl = 45^(^1 ~ V 2 + 2J)(Vi -V 2 - 2 J), (12a) 
h c2 = ^{V l + V 2 + 2J)^^_^ r (12b) 
ft c3 = 4S(Vi + V 2 + 2 J). (12c) 

In the figures, the quantities on the axes are also scaled by 
J c and ho = /i C 3|./=o to compare the results in the same 
scale. In addition to the CS and CSS phases, other solid 
(with p = 1/4, and 3/4) and supersolid (SS2a) phases are 
formed due to the competition of the NN and NNN re- 
pulsions^^ These phases have the sublattice structures 
depicted in Fig. [3] 

As seen in Eqs. (Ti"2")) , the CSS phase can emerge as long 
as the NNN interaction is finite, and the window h c i < 
\h\ < h c2 gets wider as V 2 increases. Thus, it appears 
that we just have to prepare the system with a stronger 
NNN interaction in order to obtain the stable CSS phase 
in a wider range of the parameters. However, when the 
value of V 2 /Vi is large, the striped solid order shown in 
Fig. [UlI) is more favored than the checkerboard. The 
general expression of the MF energy for striped phases is 
given by 

E { st) /M = -JS 2 (sm0 Rl +sm8 R2 ) 2 

+V 1 S 2 (cos0 Rl + cos6» R , 2 ) 2 /2 
+2V 2 S 2 cos6» Rl cos6» R2 
-hS(cos9 Rl +cos0 R2 )/2, (13) 

where 9 Rl and (9r 2 are the canting angles of the pseu- 
dospins on even and odd rows [see Fig. [TJII)] . For ex- 
ample, let us consider the solid orders emerging at the 
half-filling (p = 1/2). Putting cos6>a = — cosfe = 1 in 
Eq. pT|) and cos 8 Rl = — cos 6*r 2 — 1 in Eq. (jTB")) , we ob- 
tain the MF energies of the CS and striped solid states: 

E { ch) /M = -2S 2 (^i - V 2 ) (CS), (14a) 
E^/M = -2S 2 V 2 (striped solid). (14b) 

From the comparison, one finds that the striped solid 
state has lower energy than the CS state when V 2 /Vi > 
1/2. Also for the supersolid phase, the striped one takes 
the place of the CSS phase in this regime4^ Because of 
the transitions to the striped phases, we cannot extend 
the CSS region by exceeding the limit of V 2 /Vi = 1/2. 
Moreover, as shown in Fig.[5Jb), when the value of V 2 /Vi 
approaches the boundary to the stripe regime, the SS2a 
phase is extended toward the large J jV\ region due to 
the competition of the two density-wave orders. This 
competition also causes strong quantum fluctuations that 
destabilize the CSS states, as will be discussed in Scc. lIVI 
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FIG. 2: Ground-state MF phase diagrams of the V1-V2 model in the (J, h) plane for (a) V2/V1 = 0.2 (after Pich et al., Ref. |43T) 
and (b) V2/V1 = 0.4. Second- and first-order phase transitions are indicated by the thin and thick solid lines, respectively. The 
dashed vertical lines mark the location of J/J c = 0.85. The axes are scaled in two different ways: by Vi for both axes and by 
J c and ho, respectively. 



B. The MF results for the V d i P model 

Next, let us move onto the Vdi p model. In Ref. [I?], we 
have applied the MF theory to this model, and exam- 
ined the stability of superflow in the CSS state. Here, 
we present more detailed information on the MF ground 
states, and discuss the comparison with the results for 
the V1-V2 model. 

The MF energy per site of the Vdi p model for the 
checkerboard pattern can be written as 

E { ch) /M = ~AJS 2 sin A sin 8 B + 2Vf s S 2 cos 9 A cos 8 B 

+y 2 cff s ,2 (cos 2 e» A + cos 2 e B ) 

-hS(cos6 A +cos6» B )/2. (15) 
Here, (V^) is just the summation of the long-range 



(l)SS2a (E) P = m Solid 




FIG. 3: Schematic pictures of the (I) SS2a and (II) p = 1/4 
(p = 3/4) solid phases emerging in the phase diagram of the 
V1-V2 model. There are two possibilities (a and b) for the 
structure of the p = 1/4 (p = 3/4) solid phase, which are 
energetically degenerate for the V1-V2 model within the MF 
theory. 



interactions between the pseudospins on the same (dif- 
ferent) sublattice sites: 

V ? S = ^E^Ate=l-460V, (16a) 
V i S = ^E y ^ A = 0.7985V. (16b) 

Ia 

The index j A (j B ) means the jth site on sublattice A 
(B). Only by replacing V\ and V2 with V° s and 
in Eq. (jlip . we can immediately obtain the expression 
of Eq. (fT5"T) . This means that the MF properties of the 
checkerboard phases of the Vdi p model can be described 
exactly by the V1-V2 model with the effective NN and 
NNN interactions V* s and V 2 eS . For example, the phase 
boundaries between the CS, CSS, SF, and MI phases are 
obtained by replacing V\ and V2 in Eqs. (fT2)l with Vf s 
and V 2 cS ■ Moreover, the tip of the CS lobe is given by 
Jc = (Vf s - V 2 cS )/2. It is worth noting that the MF 
energy of Eq. (fT5|) is valid not only for the Vdi p model, 
but generally for systems with checkerboard sublattice 
structure regardless of the form of Vji . 

The resulting MF phase diagram shown in Fig. @] has a 
similar structure to that of the V1-V2 model in Fig. [21 es- 
pecially for the region of J/V > 0.2. However, many ad- 
ditional phases emerge for the region of smaller J/V due 
to the long-ranged character of the dipole-dipole interac- 
tion. Within our analysis (see Appendix IA")) . we found 
the supersolid, named SS2b, and the solid phases with 
p = 1/3 and 2/3 in addition to the phases appearing in 
the V\-V 2 model. Unlike the V\-V 2 model, the two possi- 
ble structures of the p = 1/4 (3/4) solid state shown in 
FigGHn) can be distinguished even within the MF theory; 
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FIG. 4: The same as in Fig.[2]for the V dip model. The p = 1/4 
(3/4) solid state here has the b-type symmetry in Fig. [3^11) . 
The lower panels are the sketches of the SS2b and p — 1/3 
(p = 2/3) solid states. Many other phases with more complex 
structure can emerge in the small J/V (shaded) region. The 
dashed vertical lines mark the location of J/J c = 0.6 and 
0.85. 



the b-type structure has lower energy. The emergence of 
these solid phases is consistent with the QMC results*^ 
Although many other phases can emerge for smaller J/V, 
we do not extend the calculations to more complex sub- 
lattice structures, since our main focus is the stability of 
the CSS phase. 

It should be noted that the ratio of the effective NNN 
interaction strength to the NN one is fixed in the Vdi p 
model a,sV 2 cfi /V° s sa 0.55. This value obviously exceeds 
the limit V2/V1 = 1/2, above which the striped phases 
emerge in place of the checkerboard ones in the case of the 
V1-V2 model. Nevertheless, we have to keep in mind that 
the effective interactions Vff are made by the summation 
of the long-range interactions between various pairs with 
different distances. Therefore, the limit predicted for the 
V1-V2 model cannot be directly applied to the Vdip model. 

In the Vdip model, the MF energy per site for the 



striped phases is written as 

E { Q st) /M = -JS 2 (sin Rl + sin R , 2 ) 2 

+VfS 2 (cos Rl + cos R2 f /2 

+2VfS 2 cos0 Rl cos0 Ra 

-hS(cos 6 Rl + cos 6> R , 2 ) /2. (17) 

This expression is formally equivalent to that of the V1-V2 
model [Eq. (fT3"|) ] with the effective interactions 

s = 2.025V, (18a) 

^Rl 

^° ff = jEfe 2 fe 1 -^ Rl / Rl )= 0.2339F.(18b) 

However, the effective NN and NNN interactions have 
different values for the checkerboard {V°% ) and striped 

(V12) phases [compare Eqs. (fll))) and (fT5))]. Hence, the 
large value of V r 2 eff /V 1 eff in the checkerboard phases does 
not mean that the striped phases are energetically pre- 
ferred, and the checkerboard order is always favored over 
the striped one in the Vdi p model. As an example, we 
show the comparison of the MF energies of the CS and 
striped solid states: 

E { ch) /M = -2S 2 {Vf -Vf) 

= -0.3307V" (CS), (19a) 
E^ Bt) /M = -2S 2 Vf 

= -0.1169V (striped solid). (19b) 

According to Eqs. ([12]). the CSS region, h c i < \h\ < 
h C 2, gets wider for a larger value of V2/V1 {V£ s /V° s ). In 
the Vdip model, the ratio V£ s /V* s « 0.55 is larger than 
V2/V1 of the V1-V2 model with the checkerboard order, 
the CSS region is also larger in the MF level. This is 
one of the two main reasons why the CSS phase is stable 
in the case of the dipolc-dipolc interactions. Comparing 
the width in units of ho, for example, at J/J c = 0.85 
in Figs. [2ja), [2jb) , and HI we indeed see that the Vu p 
model has a wider region of the CSS phase than the V1-V2 
model. Moreover, despite the large value of V 2 eS /Vi S , the 
SS2a region in Fig. 0] is relatively suppressed compared 
with that in Fig. [2jb). This means that the direct com- 
petition of the checkerboard and striped density- wave or- 
ders is much weaker than the case of the V1-V2 model. 
The suppression of the competition can be also seen in 
the excitation spectra, which will be discussed in the next 
section. 



IV. LINEAR SPIN- WAVE ANALYSIS 

In this section, we discuss the strength of quantum fluc- 
tuations around the MF ground states within the linear 
spin- wave (LSW) theory. 42 ' 43 ' 46 First, we perform local 
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rotations of the spin reference frame in Eq. ((4J , so that 
the new spin quantization axis is oriented along the di- 
rection of the classical pseudospin vector: 

f Sf\ ( cos9 3 sinflj \ [SJ\ 
S? = 1 \[sy . (20) 

V J V ~ j cos ^ / V / 

Furthermore, we introduce new bosonic variables via the 
Holstcin-Primakoff transformation, 

5? = S-b]bj, (21a) 
= 1 ^25-6]6A +^V25-^) , (21b) 




, (21c) 



to describe quantum fluctuations around the classical 
spin angles. Within the LSW approximation, we keep 
the terms up to the second order in the boson operators: 

-Hspin ~ E ± H 2 -, (22) 

where Eq is identical to the MF energy given by Eq. ((5J . 
The linear term in boson operators disappears by substi- 
tuting the MF solutions into 9j. Diagonalizing H 2 , we 
calculate the LSW excitation spectra w(q) and the num- 
ber of "spin waves" (bjbj) to estimate the strength of the 
quantum fluctuations (see Appendix [B] for details). By 
calculating the number of spin waves, one can roughly es- 
timate the strength of quantum fluctuations around the 
MF solutions obtained in Sec. IIIII In the spin language, 
the value of (bjbj) corresponds to the spin reduction from 
its classical value S due to the zero-point fluctuations. 

We will show the results of the excitation spectra 
in Sec. IIV Al and of the the number of spin waves in 
Sec. II V B] For the V1-V2 model, the LSW excitation spec- 
tra of SF, CS, and CSS states have already been discussed 
in detail in Ref. and it was confirmed that the soften- 
ing of roton excitations causes the phase transition from 
the SF to CSS state. As for the Vdi P model, although 
we used in Ref. Hil the LSW theory to discuss the criti- 
cal velocity of flowing CSS states, detailed results of the 
spectra have not been presented yet. Moreover, to date, 
no studies estimating the strength of quantum fluctua- 
tions from the values of (bjbj) have been demonstrated 
for the comparison of the two models. 

A. The excitation spectra 

As mentioned above, the LSW excitation spectra for 
the V1-V2 model have already been analyzed in Ref. HH. 
Hence, we show here only the case of the dipole-dipole in- 
teraction given in Eq. <j3j> . In the calculations, we have to 
take an infinite summation in Eq. (|B12j) due to the long- 
range nature. To avoid the practical difficulty, we intro- 
duce a cutoff distance on the dipole-dipole interaction as 



(a) SF (b) CSS 




FIG. 5: Excitation spectra w(q) of the Vdi p model in the (a) 
SF (at h/ho = ±1), (b) CSS (at h/ho = ±0.58), and (c) SS2a 
(at h/ho = ±0.5) phases for J/J c = 0.6 . 



Vji = for \rj — r;| > I6d only in this section (namely, 
in Sees. IIV Al and HV B|) . For this truncated dipole-dipole 
interaction, the values of the effective NN and NNN in- 
teractions are = 1.410V and V 2 cB = 0.7495V and 
the ratio is V 2 eS /Vf s « 0.53. 

Solving Eq. (|B14|I . we plot in Figs. E£a)-G^c) the ex- 
citation spectra w(q) for the SF, CSS, and SS2a phases 
along the line J/J c = 0.6 marked in Fig. [4] The excita- 
tion spectra have a Nambu-Goldstone mode reflecting the 
spontaneous breaking of the U(l) symmetry. The spec- 
trum of the CSS phase consists of two branches due to the 
two-sublattice structure, and the lower branch has gap- 
less, linear dispersions around q = (0,0) and (n/d,ir/d), 
which is the ordering vector of the checkerboard phases. 
In the SS2a phase, in which the checkerboard and stripe 
orders coexist, the lowest mode is gaplcss at q = (w/d, 0) 
in addition to at q = (0,0) and (n/d,n/d). 

Figures |U[a)-[fJd) show the excitation spectra at the 
phase transitions between the different phases. When 
one approaches the CSS phase from the SF region, a 
roton-like minimum at q = (jr/d,ir/d) develops, and it 
touches zero at the boundary h = ±h c2 as shown in 
Fig. Eta), causing the second-order phase transition to 
the CSS state. In a similar way, a roton-like mode at 
q = (ir/d, 0) causes the second-order transition from CSS 
to SS2a [see, Fig.[SJb)]. In contrast, as shown in Fig.^c), 
such a signal does not appear in the spectra at the first- 
order phase transitions. The above-mentioned properties 
of the excitations qualitatively agree with the case of the 
V1-V2 model. 
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(a) from SF to CSS 
(second-order) 



(b) from CSS to SS2a 
(second-order) 




(0,0) (1,0) (1,1) 

qd/n 



(0,0) (0,0) (1,0) (1,1) 

qd/n 



(0,0) 



(c) between CSS and SS2a 
(first-order) 

1 SS2a 




(0,0) (1,0) (1,1) 

qd/n 



(0,0) (0,0) (1,0) (1,1) 

qd/n 



(0,0) 



FIG. 6: Excitation spectra oj(q) of the Vdi p model at the 
second-order transition (a) from the SF to CSS phase (at 
h — ±/i C 2) and (b) from the CSS to SS2a phase (at h/ho = 
±0.3819) for J/J c = 0.6. (c) The same as in panels (a) and 
(b) at the first-order transition between the CSS and SS2a 
phases (at h/ho = ±0.5348). The CSS and SS2a states are 
energetically degenerate at this point. 



B. The number of spin waves 

For all the cases of Figs. (UJa), HJb), and[U the system 
exhibits the phase transition from SF to CSS, and then it 
reaches the CS phase if the value of h/ho increases from 
a negative value to zero along the line of J/J c = 0.85. 
We will plot the number of spin waves (bjbj) along this 
line as a function of the filling factor p. Within the MF 
analysis, the filling factor p shows a linear increase with 
the chemical potential h both in the SF and CSS phases, 
and the slope of the line, which is proportional to the 
compressibility, is always larger in the CSS phase than 
in the SF phase: 



r 1/2 

P= 1/2 



h/2h c3 (SF) 
(h T h cl )/l6V 2 S (CSS) 



(23) 



where the upper (lower) signs are for positive (negative) 
values of h. Recall that V2 has to be replaced with V 2 cS 
for the Vdip model. It should be noted that the critical 
filling factor p c at the SF-CSS transition point, which is 
obtained by substituting h = ±/i C 2 into Eq. (|23[) . can be 
expressed as a function only of J/J c : 



= I±I / — J l Jc 

Pc 2 2yi±j/j c 



(24) 



It takes p c = 0.3576 for J/J c = 0.85 in the low-density 
(negative h) side. 

Figures [7fa)-|7fc) show the results for the V1-V2 model 
with V2/V1 = 0.2 and 0.4, and for the Vdi p model. In all 

the cases, we can see that (bjbj) has a peak at the SF-CSS 
phase transition, which means that quantum fluctuations 
are particularly strong at the phase boundary. The num- 
ber of spin waves for V2/V1 = 0.4 in Fig. [TJb) is much 
larger than the case of V 2 /Vi = 0.2 in Fig. Efa). This 
is attributed to the strong competition of the NN and 
NNN interactions. In fact, as shown in the right panel of 
Fig. Elk), the excitation spectrum w q exhibits a remark- 
able drop at q = (ir/d,0), which indicates the existence 
of strong striped density-wave fluctuations. The maxi- 
mum value of (bjbj) for V2/V1 = 0.4 reaches about 77 
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(b) F 2 /Fi = 0.4 CS 
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0.4 0.5 
'ft 
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(c) dipole, Kf ff /Kf ff =0.53 

(up to distances of 16J) 



to 0.2 



SF 


CSS 







0.0 0.1 0.2 0.3 ^0.4 0.5 

P v c 




(0,0) (1,0) (1,1) (0,0) 

qd/n 



FIG. 7: The number of spin waves (bj a bj a ) as a function of the 
filling factor p for J/J c = 0.85 and the excitation spectra w q 
at the SF-CSS phase transition point p = p c . The panels (a), 
(b), and (c) show the results for the V1-V2 model with V2/V1 = 
0.2 and 0.4, and for the Vdip model with the truncated dipole- 
dipole interaction, respectively. In the CSS phase, (bj a bj a ) 
takes two different values on each sublattice. 
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percent of the classical value of the spin length S = 1/2, 
which means that the predictions for the ground states 
within the MF theory are unreliable. Actually we will 
show in Sec. [V] that the CSS phase predicted by the MF 
theory almost completely disappears due to the strong 
quantum fluctuations. 

On the other hand, the excitation spectrum w q for the 
Vdip model in the right panel of Fig. [7|c) does not ex- 
hibit a significant drop at q = (ir/d, 0) unlike the case 
of V 2 /Vi = 0.4. This comes from the fact that V£ s is 
not just the NNN interaction but the summation of vari- 
ous long-range interactions that weakens the competition 
with the stripe order. Therefore, as shown in Fig. [TJc), 
the quantum fluctuations in the Vdip model are relatively 
weak for its large value of V£ s /Vj* . This is the second 
reason for the stability of the CSS state in the Vdip model. 
The CSS region predicted by the MF theory is signifi- 
cantly reduced by the quantum fluctuations but still re- 
mains sufficiently large (see Ref. [TH and Sec. \V\ of this 
paper). 



V. LARGE-SIZE CLUSTER MEAN-FIELD 
METHOD AND SCALING ANALYSIS 

We drew the ground-state phase diagram of the two 
models within the MF theory in Figs, E^a), C^b), andg] 
However, according to the previous section, the fluctua- 
tions around the classical ground states are too large to 
completely ignore in any case. In this section, considering 
the MF results obtained in Sec. IIHI as a starting point, 
we discuss how the quantum fluctuations change the fea- 
tures of the ground-state phase diagrams by employing 
a large-size CMF method^ We perform the calculations 
based on rectangular-shaped clusters, and then extrap- 
olate the results with respect to the cluster size. The 
obtained results will be compared with the QMC data 
in Ref. [H for the Vdi p model. As for the V1-V2 model, 
although the authors of Ref. H3 concluded that the CSS 
state is thermodynamically unstable from the QMC cal- 
culations for V\ = 3 J, the entire phase diagram including 
the effects of quantum fluctuation has not been produced 
yet. We will confirm, in the (J/Vi, ft./Vi)-plane, that the 
CSS phase is almost completely destroyed by the strong 
quantum fluctuations in the V1-V2 model. Hereafter, we 
will use the full (untruncatcd) dipolc-dipole interaction 
again in the Vdip model. 



A. The CMF method 

First, we describe the details of our CMF approach. 2 - 
The standard MF theory approximates the system by 
single-site problems in effective fields. A natural exten- 
sion of the single-site approximation is the use of "clus- 
ters" of multiple sites as an approximate system4Zr— For 
example, the Bcthc-Pcicrls- Weiss (BPW) methods- 
employs a cluster consisting of one central site and its 



(I) 4x4 cluster 




A OB 



FIG. 8: Clusters of (I) 4 x 4 sites and (II) 3x3 sites embedded 
into the checkerboard sublattice pattern. 



directly connected sites, e.g., a cluster of (1 + 6 + 6) sites 
for a triangular-lattice system with NN and NNN inter- 
actions^ Treating exactly the interactions within the 
cluster, one can partially take into account the effects of 
correlations between particles (or spins). However, the 
BPW method and its extension a 51 ' 52 cannot be applied 
to an infinite-range interaction model like the Vdi p model 
since all sites are "directly connected" by the long-range 
interactions. 

Oguchi's method^ is another simple way to extend 
the MF theory to clusters. In Ref. H3, Oguchi stud- 
ied fcrromagnctism and antiferromagnctism of the low- 
dimensional Heisenberg model by using a cluster of up 
to three spins to include the short-range correlations be- 
tween the spins. Since the influence from the spins out- 
side of the cluster is also included as effective internal 
fields, we can treat even a system with infinite-range in- 
teractions. However, as Oguchi himself pointed out, the 
cluster of two or three sites is too small to sufficiently take 
into account the effects of the correlations (or quantum 
fluctuations) . 

Our CMF approach 2 ^ is an extension of Oguchi's 
method to larger-size clusters and to multiple-sublatticc 
problems. Although we use here the pseudospin form 
of the Hamiltonian, Eq. to explain the procedure of 
our method, the same manner can be applied straight- 
forwardly to hardcore bosons and even to softcore-boson 
models. First, we assume a sublattice structure expected 
to emerge in the parameter range. Then, we embed a 
cluster of Ac sites into the background sublattice struc- 
ture. Figures [Sjl) and [8)^11) show, as examples, the cases 
of Ac =4x4 and Ac = 3x3, which we refer to as 
CMF-4 x 4 and CMF-3 x 3 under the assumption of the 
checkerboard sublattice structure. As for the case of the 
CMF-3 x 3, we have two inequivalent choices for embed- 
ding the cluster. We have to deal with both of them 
equally as in the BPW method for multiple-sublatticc 
systems^. Now, instead of treating the many-body prob- 
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lem in the whole system given by Eq. ((4]) , we consider the 
effective cluster Hamiltonian Hq written as the following 
general form: 

H C = -23 (?f S i + Spf) V ^f S l 

{j,i)ec j.iec 

-E( /i + ^ cff )^-E^ 0ff ^> (25) 

in which the interactions within the cluster are treated 
exactly, while the interactions between the spins in the 
cluster and the rest of the system arc approximately in- 
cluded via the effective fields 
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and m^'' 

the CMF method, which act as the mean fields from the 
spins in C . Here, we chose again the xz plane as the 
plane in which the spins lie; i.e., (Sj) = 0. If we have 
two or more possibilities of choosing the cluster like in 
the CMF-3 x 3, we should consider all the corresponding 
cluster Hamiltonians like Hc 1 , Hc 2 > " ' ■ 

Note that in our CMF method, we consider the Nq- 
site problem in the cluster just as a reference system to 
estimate the values of the mean fields rn^ x , which depend 
only on the background sublattice index of the site; i.e., 
mj" ee mf' x . For example, the effective fields acting on 
the top-left site "1" in the 4x4 cluster of Fig. Ejl) can 
be written as the following explicit forms: 



0.15 0.20 0.25 0.30 

J IV 



z,eS 



-214m! - 3V 2 m% 



A- 



4 Jmg 



(27) 



for the V1-V2 model and 



7 z.ctf 




A3m^ 



for the Vdip model (see Appendix [C] for more details) . 
The values of the mean fields m z a and m x a are calculated 
sclf-consistcntly as the expectation values of the pseu- 
dospins inside the cluster as follows: 



m a X = (S^ x )cmf 



M c N a ^ 



Tr [s z 3 - a x e-' 3H ^ 
Tr (e-^- ) 



(29) 



FIG. 9: Phase boundaries between the CS, CSS, and SF 
phases by the CMF-4 x 4 calculations for the V1-V2 model 
with (a) V2/V1 = 0.2 and (b) V2/V1 = 0.4, and for (c) the V d ip 
model. The thin and thick solid lines indicate the second- 
order and first-order transitions. For comparison, the MF 
results, Eq. (|12p . are shown by the dashed lines. 



where j3 = 1/T (we take T — > in this paper), Mq is 
the number of the possible choices of the cluster, and 
N a ee ^ n N a ^ n is the summation of the number of 
sites belonging to sublattice a in cluster C n . For ex- 
ample, Mc = 2 in the CMF-3 x 3 for checkerboard 
phases in Fig. [H(II) and we have (Na,i,Nb,i) = (5,4) 
and (N\ t 2, -^8,2) = (4, 5) for the two clusters C\ and C2, 
respectively, which leads to Na = Nb = 9. 

This method reduces to the conventional MF (namely, 
Weiss's molecular- field) theory for Nc = 1, and becomes 
exact in the limit Nc — > 00. We have to diagonalizc 
the cluster Hamiltonian to take the trace on the right- 
hand side of the self-consistent equation, Eq. (f2"!)| . Thus 
the practical limit of the cluster size Nc is determined by 
the largest number of sites which can be treated by exact 
diagonalization techniques. It should be noted, however, 
that some of the symmetries of the original Hamiltonian, 
Eq. (j4]), are broken in the effective cluster Hamiltonian 
due to the existence of mean fields. 



B. The CMF results for the 4x4 cluster 

First, we show the results of the CMF-4 x 4 method for 
the two models and compare them with the MF results. 
We focus on the checkerboard phases and discuss the 
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influence of quantum fluctuations on the locations of the 
phase boundaries between the CS, CSS, and SF phases, 
namely h — ±/i c i and ±/i C 2 m Eqs. (fl2|) . Note that there 
are no quantum fluctuations at the boundary between the 
SF and MI phases, h = ±/i C 3, and thus the expression in 
Eq. (fT2b) does not change for any Nq- Other phases 
with more complex symmetries, such as three-sublattice 
and four-sublattice phases, are affected more strongly by 
the quantum fluctuations, and the locations are shifted 
towards the region of much smaller values of J / J c than 
those of the MF theory . 18 ! 57 Thus, we leave them out 
of the scope of the rest of this paper, treating only a 
relatively large-J/ J c region. 

For the checkerboard phases, we can obtain the SF 
order parameter \^\ = \m\ + m%\/2, the CS order pa- 



rameter \pq\ 



/2, and the filling factor p 



1/2 + (m z A + mg)/2 by solving Eq. (j29j) self-consistcntly. 
The second-order transition boundary from CSS to CS 
is determined by the point at which |>I'| vanishes (or at 
which the value of p reaches 1/2). On the other hand, the 
CS order parameter \pq\ vanishes (namely, m z A = m^) 
at the second-order transition from CSS to SF. In ad- 
dition to these second-order (continuous) transitions, we 
find that first-order (discontinuous) transitions between 
CS and SF and between CSS and SF can also appear 
due to the effects of the quantum fluctuations. In the 
CMF formalism, we cannot directly calculate the value 
of free energy of the system. Instead, we use the Maxwell 
construction in the (J, x) plane to determine the first- 
order phase boundaries. The quantity x is defined by 

Figures [9fa)-[9jc) show the phase boundaries between 
the CS, CSS, and SF phases obtained by the CMF-4 x 4 
method with the corresponding MF results. In all the 
cases, the regions of the CS and CSS phases shrink con- 
siderably because of the quantum fluctuations. The re- 
ductions of the values of J c from the MF values are 5.5 
percent for V2/V1 = 0.2, 16.5 percent for V2/V1 = 0.4, 
and 11.1 percent for the Vdi p model within the CMF-4 x 4 
level. This fact indicates that the quantum fluctuation 
of the V1-V2 model is stronger for a larger value of V2/V1, 
and that of the Vdi p model is small relative to the large 
value of V 2 eS /Vi S , which is consistent with the spin- wave 
prediction in Sec. IIVI In Fig. 1101 we compare the reduc- 
tions of the value of \h\ at the CSS-SF boundary, h c2 , 
from the MF value given in Eq. (|12b ). This result also 
confirms the statement made in Sec. IIVI 

Other than the shrinking of the CS and SS phases, we 
can see a qualitative difference between the CMF and MF 
results. As clearly seen in the enlarged views, Figs. fTlT a) 
and lllf b). the direct first-order (discontinuous) transi- 
tion from the CS to the SF phase emerges in the CMF 
result, and the transition between the CSS and SF phase 
also becomes discontinuous near the triple point of the 
CS, CSS, and SF phases. This result is attributed to the 
fact that the classical degeneracy of the CS, CSS, and SF 
states at J = J c and h — is lifted by taking into ac- 
count the quantum fluctuations. However, the first-order 
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FIG. 10: The ratio of h c2 obtained by the CMF-4 x 4 method 
to the MF value at fixed J/ J c ■ 
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FIG. 11: Enlarged views of the region around the tip of the 
lobe in Fig.H;b) [panel (a)] and Fig.^c) [panel (b)]. The first- 
order transitions from the CS to SF phase and from the CSS 
to SF phase are found. Also in the case of V2/V1 = 0.2, we 
find a narrow but finite region where the first-order transitions 
occur. 



transitions occur in a narrow range of parameters and the 
discontinuity of the order parameter is very small, and 
thus the first-order nature has not been reported previ- 
ously in the QMC works.^ 7 - 



C. The cluster-size scaling 

We perform the infinite-size extrapolation, Nq — > 00, 
of the CMF results with different-size clusters. We use 
a series of rectangular-shaped clusters of Nq = 1x2, 
2 x 2, 2 x 3, 3 x 4, and 4x4, which are shown in TablelU 
The clusters with odd numbers of sites, e.g., Nq = 3x3, 
are not treated here, because they may belong to a dif- 
ferent scaling series from that of the clusters with even 
numbers of sites. To perform the infinite-size extrapo- 
lation, we introduce the scaling parameter A defined by 
Nc^xz/2 > wmc h varies from to 1. Here, N-q is the number 
of bonds within the cluster and the denominator means 
the number of bonds of the original lattice per Nq sites. 
The parameter A provides an indication of how much the 
correlation effects between the particles are taken into ac- 
count by using the cluster. The value of A for each cluster 
is listed in Table fl] Note that the MF (JV C = 1) and exact 
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(N c = oo) results correspond to A = and A = 1, respec- 
tively. The accuracy of the scaling procedure (A — > 1) is 
discussed in Appendix iDl 

Now we perform the scaling analysis to the three cases, 
V 2 /Vx = 0.2, Vv/Vx = 0.4, and the V dip model. We first 
consider the change in the location of the tip of the CS 
lobe (J = J c ) with increasing Nc- In all the cases, the 
value of J c /Vx (or J c /V for the Vdi P model) systemati- 
cally decreases with the cluster size Nc, and the linear 
fits of the data for the three largest clusters (Nc = 2x3, 
3x4, and 4x4) are fairy good as shown in Fig. IT27 a). 
The scaled values ofJJVx (J c /V) are 0.3701, 0.2318, and 
0.2817 for V 2 /Vx = 0.2, V 2 /Vx = 0.4, and the V d ip model, 
respectively. Next, we move on to the scalings of the 
phase boundaries between the CS, CSS, and SF phases. 
The extrapolations are carried out on the value of h/ho 
at each transition at fixed J/ J c , in which J c is the value 
at each cluster size. We show examples of linear fittings 
of the CMF data for the phase boundaries between the 
CS and CSS phases (h = ±/i c i) and between the CSS 
and SF phases (h = ±h c2 ) in Fig. IT2Tb). We can see 
that in the cases of the Vx-V 2 model, the two lines of 
the CS-CSS and CSS-SF transitions approach each other 
very closely in the limit A — > 1. Especially, the lines for 
V 2 /Vx = 0.2 intersect before reaching A = 1, which means 
that the transitions are replaced by the direct first-order 
transition between the CS and SF phases. On the other 
hand, the CSS region (h cl < \h\ < h c2 
remains sufficiently large. 

Performing the same scaling analyses on the CS-CSS 
and CSS-SF (or CS-SF) transition boundaries for other 
values of J/J c , we draw the expected phase diagrams for 
the limit A — >• 1 in Figs. IT5ra)-IT37c) . Here, the quantities 
on the axes are rescalcd in units of V\ (or V) by using 
the scaled value of J c /V\ (J c /V), obtained in Fig.JT^a), 
and ho/Vi (h n /V) for each case. We can see that the 
width of the CSS phase almost vanishes in the two cases 
of the Vi-V 2 model. Within the accuracy of the scaling 
procedure, it is difficult to provide a final conclusion on 
whether a very small region of the CSS phase can survive 



TABLE I: A series of clusters used in our CMF calculations. 
The values of Nb and A are also listed. 



of the Vdip model 
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FIG. 12: (a) Cluster-size scalings of the CMF data for the 
value of Jc/Vi (J c /V). The lines are the linear fits of the three 
points, A c = 2 x 3, 3 x 4, and 4x4, for V2/V1 = 0.2 (triangle, 
dashed line), V3/V1 — 0.4 (square, dash-dotted line), and the 
Vdip model (circle, solid line), (b) Cluster-size scalings of the 
CMF data for the phase boundaries between the CS and CSS 
phases, h = ±/i c i (open symbols), and between the CSS and 
SF phases, h — ±/i C 2 (closed symbols), at J/J c = 0.7. The 
triangle, square, and circle symbols correspond again to the 
data for V2/V1 = 0.2, V2/V1 = 0.4, and the V dip model. 



or completely disappear. However, even if the CSS region 
can survive, it should be too narrow to detect, and this 
result does not contradict the conclusion of Ref. H3- Fig- 
urc fTST c) shows the scaled CMF result for the Vdip model, 
which is in surprisingly good agreement with the QMC 
datapi£ Unlike the two cases of the V1-V2 model, we can 
see that the CSS phase remains stable in a considerably 
large region of parameters. 

Having obtained the above CMF results, we now sum- 
marize the difference of the V1-V2 and Vdip models to- 
gether with the knowledge gained from the MF and LSW 
analyses in Sec. Mil and Sec. IIV1 As for the V1-V2 model, 
we have the following dilemma: the NNN interaction V 2 is 
required to be large in order to obtain a large region of the 
CSS phase, according to the MF prediction in Eq. (fT2"l) ; 
however, the LSW analysis showed that the larger the 
value of V2/V1 is, the stronger quantum fluctuations are 
due to the competition between the checkerboard and 
stripe density-wave orders. Because of this dilemma, the 
CSS state cannot be stabilized in a sufficiently large re- 
gion of the phase diagram for both cases of small and 
large values of V2/V1, as shown in Figs. H3fa) andfTffib). 
On the other hand, the MF phase diagram of the Vdi p 
model contains a rather large region of the CSS phase 
thanks to the large value of V 2 eS /Vf s and at the same 
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FIG. 13: The results of the scaling analyses to the phase 
boundaries between the CS, CSS, and SF phases for the Vi- 
V 2 model with (a) V2/V1 = 0.2 and (b) V2/V1 = 0.4, and 
for (c) the Vdip model. First-order transitions between the 
CS and SF phases are expected to occur at the thick solid 
lines and in the region between the two dotted lines, which 
correspond to the "metastability limits" of the CS and SF 
phases. For comparison, the MF results, Eq. (JX3J) , and the 
QMC data (Ref. Ila ) are shown by the dashed lines and by 
the squares and circles, respectively. 



time, the long-range nature of the dipole-dipole interac- 
tion suppresses the quantum fluctuations around the MF 
ground state because the competition with other solid or- 
ders is weaker. For these reasons, the CSS phase can sur- 
vive in the Vdi p model, as shown in Fig.ljUTc). even after 
taking into account the effect of quantum fluctuations. 

A similar discussion is applicable to the difference be- 
tween triangular and kagome lattices in the hardcore 
Bose-Hubbard model with nearest-neighbor interaction. 
The MF properties of the two systems are identical ex- 



cept for the scale of the chemical potential, and both the 
systems have a large region of supersolid phase in the 
MF phase diagram^ . However, previous QMC studies 
have shown that while a stable supersolid phase exists in 
the triangular-lattice system^ it has not been found in 
the kagome latticed This is a similar situation to the 
difference between the V1-V2 and Vdip models. In this 
case, although the competition of different solid orders 
does not make a large difference between the two sys- 
tems, it is known that the quantum fluctuations in the 
kagome lattice are much stronger than those in the trian- 
gular lattice, reflecting, e.g., the lower coordination num- 
ber^ Because of the strong quantum fluctuations, the 
supersolid states in the kagome lattice are more strongly 
destabilized and cannot survive in the QMC calculations. 
Therefore, we can say that for the emergence of stable lat- 
tice supersolid states, in general, it is necessary to satisfy 
the following two (qualitative and quantitative) condi- 
tions: A certain long-range interaction Vji is required for 
creating solid orders, and the quantum fluctuation should 
be weak enough so as not to destabilize the supersolid 
states into phase separation. 



VI. SUMMARY 



In conclusion, we have investigated the ground-state 
phase diagrams of the hardcore Bose-Hubbard model 
with square lattice structure and long-range interactions. 
One of our main focuses is placed in understanding 
the role of long-range interactions in the emergence of 
checkerboard supersolid (CSS) states, through the com- 
parison of the models with nearest-neighbor and next- 
nearest-neighbor interactions (the V1-V2 model) and with 
the dipole-dipole interaction proportional to 1/r 3 (the 
Vdip model). Specifically, we discussed the reasons why 
the CSS states can be stable only in the case of the Vdip 
model, and clarified the origin of the qualitative differ- 
ence between the two systems. We first showed the clas- 
sical (mean-field) properties of the systems, and then 
discussed the strength of quantum fluctuations around 
them in terms of the linear spin-wave theory. Moreover, 
we also applied the cluster mean-field (CMF) method 
and its cluster-size scaling to take into account the ef- 
fects of quantum fluctuations in a self-consistent way. 
We confirmed quantitative accuracy of our CMF scaling 
procedure^ by making a comparison with the quantum 
Monte Carlo data in Fig. fTSt c) and Appendix Q3] In prin- 
ciple, this approach can be also applied to any other or- 
dered systems including softcore bosons and higher-spin 
systems. Especially, our CMF method may be useful 
for studying frustrated systems since it is free from the 
minus-sign problem. 
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Appendix A: three- and four-sublattice structures 

In the case of the Vdi p model, in addition to the uniform 
(SF and MI) and checkerboard (CS and CSS) phases, 
many different types of solid and supcrsolid phases can 
appear due to the long-range nature of the dipolc-dipolc 
interaction. In general, we encounter various phases with 
more complex sublattice structures when the effects of 
the dipole-dipole force are stronger (the value of J/V is 
smaller) i 18 i 29 

To obtain the results in Sec. IIII1 we restricted our MF 
analysis to the phases with up to the three- and four- 
sublattice structures shown in Figs. IT47 D and[T4TlIh fo- 
cusing on the region of relatively large values of J/V 
(but still less than 1). Each state is characterized by the 
classical pseudospin angles 0„ (n = 1,2,3,4 for the four 
sublattices; n = 1,2,3 for the three sublattices) . The 
MF energies per site for the two types (a and b) of four- 
sublattice phases are written as 

E^/AI = -JS 2 (sm6 1 + sin0 3 )(sin0 2 + sin0 4 ) 

+Vi (4a) S 2 (cos 0j + cos 3 ) (cos 2 + COS 4 )/2 

+ V 2 (4a) S 2 (cos 0i COS 3 + COS 2 COS 4 ) 

4 

+V 3 (4a) S 2 ^cos 2 n /2 

n=l 

4 

-hS^cos0 n /4, (Al) 

n=l 

I 



where Vf D > = 2 J2 h {V jah - V jlh ) /z « 0.9077V, V^°> = 



where v/ 4a) = 2j2 h V Jlh /z a 1.460V, V 2 (4a) = 



(i) Four sublattices (II) Three sublattices 



a b 




FIG. 14: Schematic pictures of (I) two types of four-sublattice 
structures (4a and 4b) and (II) a three-sublattice structure. 
The lattice sites with the same number belong to the same 
sublattice. The dashed boxes are drawn to show clearly the 
sublattice structures. 

£l,W* « 0.5162V, and V 3 (4a) = £ (l « 
0.2823V, and 



(A2) 



E h ( 2V hh + V Jlh - V nh )/z a 0.7923V, and V^' = 



£< 4b) /M 



2 sin 3 sin 4 + (sin0, + sin0 2 )(sin0 3 + sin0 4 )] /2 

2 ) (cos 3 + cos 4 )] /4 



~JS 2 [2 sin 0i sin 2 

h Vi (4b) S 2 [2 cos 0i cos 2 + 2 cos 3 cos 4 + (cos 0i + cos 
^ Ah) S 2 (cos 0i +cos0 2 )(cos0 3 + cos0 4 )/2 

4 • \ 
i cos 4 + 2J cos 2 n /4 

71=1 ' 



+ V 3 (4b) S 2 (2 COS 01 COS 2 + 2 COS 3 ( 

4 

-/iS^cos0„/4, 
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2^ Zi Vj 1 i 1 /z w 0.5584^. The above expressions are for- 
mally identical to those of the hardcore Bose-Hubbard 
model with up to the third-nearest-neighbor interactions 
(7, V2, and 7s); we have only to replace 7„ (n = 1,2, 
and 3) with the effective interactions V^ 4 ^ for the a-type 



structure or Vn for the b-type structure in the corre- 
sponding expressions of the MF energies. On the other 
hand, the MF energy for the three-sublattice structure of 
Fig. [T4Tc) is given by 



4 3) /M 



-AJS 2 {s\n9 1 sin 2 + sm6 2 smd 3 + sin6> 3 sin6>i)/3 + 2V^ 3) S' 2 (cos6>i cos6> 2 + cos6» 2 cos6» 3 + cos6» 3 cos0i)/3 



COS 6*i COS 02 + COS 62 COS 63 + COS 83 



'i)/3 + 2Vf'S*(cos 0i cos 2 + cos 2 cos 3 

3 \ 3 
COS 

01 COs2 °n / 3 - hS XI C ° S 6 ' n / 3 ' 

71 — 1 / n—1 



(A3) 



where v} 3) = £^(27^-7^)/* « 1.3177 and 7* 
2£ Il W*« 0.94177. 

To obtain the ground-state phase diagram in Fig. 21 
we carried out the minimization of the MF energies in 
the standard way. First, we minimized separately the 
MF energies given in Eqs. ([Al]). (|X2]) . and (fMJ with 
respect to the angles and then compared the three 
minimized values of the MF energies. For example, if 
such a minimization scheme leads to a solution of the 
a-type four-sublattice structure with 02 = 04, it means 
that the ground state is in the SS2a phase within the MF 
approximation. As for the other phases seen in Fig. 21 
the assumption of the b-type four-sublattice structure 
includes solutions of the SS2b and p — 1/4 solid states 
and the three sublattice structure given in Fig. IT4TII) 
includes the p = 1/3 solid state. 



K3) 



Appendix B: Details of the LSW analysis 

We present here the details of the LSW calculations. 
In the rotated frame of Eq. (pZO"]) . the pseudospin Hamil- 
tonian [Eq. (|4])] is rewritten as 



MS 



j,l (t,v=x,y,z 



(Bl) 





3 3 




where 






hi ~ 


2Jji cos 0j cos 0i — Vji sin 0j sin 


, T VV - 2 Tr 

I, 


JZZ 

hi - 


2Jji sin 0j sin 0i — Vji cos 0j cos 


TZX 

hj ~ 

hi ~ 

h x — 


= 2Jji sin 0j cos 0i + Vij cos 

V - hi ~ hi ~ '■ 
—haiadj, hj = hcosdj. 


9j sin 9i, (B2) 



After performing the HP transformation in Eq. ([2~Tj) , the 
quadratic part of the pseudospin Hamiltonian is obtained 



Ho = 




(B3) 



In order to diagonalize H 2 , the use of the equation of mo- 
tion method with the Green's functions may be more con- 
venient than the usual Bogoliubov transformation tech- 
nique for the states with complex sublattice structures. 

First, we define the retarded commutator Green's func- 
tion in a matrix form as 



«b 3 -(t) ; bj(0» = -ie{t-t'){[hMM{n) (B4) 



with 




The Green's function satisfies the following equation of 
motion: 



^<<b^);b](i')>) = 5{t-t')([hj{t),h\{t')]) 



+(([b j ,^ 2 ](t);b I t (f))).(B6) 

To solve the above equation, we perform Fourier trans- 
form of the Green's function Gjiit-t') = ((b^b]" (*'))) 
into energy space 

/OO 7 
pG^e-™ (*"*') (B7) 

and into momentum space. It is to be noted that the 
number of sites in the "unit cell" varies depending on 
the sublattice structure. For example, since the lattice 
sites are divided into two square sublattices [A and B in 
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Fig- El)] for the CS and CSS states, the Fourier transfor- 
mation into the momentum space should be performed 
on each sublattice: 

G J - < ,i, = |E G ^ e " <q ' ( '**" r, " ) . ( B8 ) 

q 

I 



where the subscripts a, j3 denote the sublattice index 
(A or B), and the sum is taken over the N/2 q- values 
in the reduced Brillouin zone. Now, we can rewrite the 
equation of motion for the Green's functions in a 4 x 4 
matrix form: 



G 



AA 



G 



AL> 



pBA r"BI 



CT Z 
CT 7 



pAA pAB 
-p§A -p§B 

q q 



G^ G 



al; 



(B9) 



The 2x2 submatriccs in Eq. (|B9| are given by 



1 

-1 



and = 



r a/3 

1 q,ll 1 q,12 
i q,12 q,H 



(a,/3 = A,B), 



(BIO) 



where 



r 



AA 
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r BB 
q,ll 

pAA 
q,12 

r AB 
q,H 

r AB 
q,12 



i? cos0 A + 45(2 J sin A sin0 B - 14 cos6» A cos6> B - V 2 cos 2 6> A ) + 25 sin 2 6 A V^\ 
H cos B + 45(2 J sin 9 A sin B - Vi cos 9 A cos B - V 2 cos 2 6> B ) + 25 sin 2 6» B V q (2) , 
25sin 2 A V^ 2) , r™ 2 = 25sin 2 B K( 2 ), 

= -4J5(cos0 A cos0 B + l)7q + 25sin(9Asin6l B ^ q (1) , 
r q^2 = -4J5(cos6» A cos6l B - l) 7q + 25sin<9 A sm9 B V^l 

I 



(Bll) 



Here, 1^ and V 2 should be replaced with V£ s and V 2 eS , 
respectively, for the Vdi p model. The Fourier factor 7 q 

( 1 2) 

is given by 7 q = [cos(q x d) + cos(q y d)]/2, and Vq ' are 
defined by 



V r 



(i) 



Ey. . e iq-(r JA -r lB ) 



y(2) = i 



Ia 



-n A ) 



(B12a) 
(B12b) 



For the V1-V2 model, these can be simply written as 
Vq (1) = V"i7 q and Vq (2) = V 2 cos(q x d) cos(q y d). 

We can obtain the spin- wave excitation spectra w(q), 
which correspond to the poles of the Green's functions, 
by solving the equation 



The number of spin waves on each sublattice (a ~A or 
B) is given by 



— lim 

2tt N ^ 
q 



die 



— 1 

x [G^ u (w + i*)-G?° u (w-W)](B15) 



where G^" u is the (1, l)-componcnt of G q Q in Eq. (|BT4)) . 

We took here the case of the checkerboard phases as 
an example. The extension to other sublattice structures 
is straightforward. For instance, the calculation for the 
SS2a phase, in which the lattice sites are divided into 
four square sublattices with lattice constant 2d, requires 
the use of a 8 x 8 matrix form instead of Eq. (|B9[) . 



dct 



w(q)l 



pAA pAB 
-p§A j^B 

q q 



(B13) 



Moreover, we can calculate (bjbj) by applying the spec- 
tral theorem to the Green's functions 



qAA qAB 

^q Wl 



ul - 



pAA pAB 
Y%A p§B 

q q 



-, -1 



<T 2 

a 7 



(B14) 



Appendix C: Details of the CMF calculations 



We present here the explicit forms of the effective fields 
hj' cS and h x -' cS in Eq. (|25p for the reader's convenience. 
Here, as examples, we show the expressions in the CMF 
calculations for checkerboard phases with Nq =3x3 and 
4x4 clusters [see Figs. IT5TI) and [15^11)]. 

In the case of Nq = 4 x 4, we have only one choice, 
Fig. [T51D. for the cluster embedded into the background 
sublattice structure. In that sense, we can say that the 
use of clusters with even numbers of sites is compatible 
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with the two-sublatticc checkerboard pattern of solid or- 
ders. The explicit forms of the effective fields for site "1" 
in Fig. 11511 ) were already shown in Eqs. (|27j) and ([28)) . 
For the other sites, we have 



(I) 4x4 cluster (n) 3*3 clusters 




FIG. 15: The same as in Figs. \8[I) and EpI) with the site 
labels (1, 2, 3, • ■ ■ ) within the clusters. 



h z ' cS = h z { cS , hl]lo = -2V x m\ - 3V 2 m%, hl' e 5 911 = -Fim| - 2V 2 m z A , h%% flil2 = -F x m A - 2V 2 m z i , 

iz.eff _ n h x,cB _ ,x,otf ,x,eS _ . j x , x,eS _ cy j x uX,cS _ „ j x ,x,cS _ V^ 1 ) 

"13, 14, 15, 16 ~ u ' "7 ~ "l i ""4,10 — U'^A) "-3,5,9,11 — ZJ '"B; "2,6,8,12 ^""'A' "13,14,15,16 ~ u 

for the V1-V2 model and 
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for the Vdip model. Using these expressions, we solved the CMF self-consistent equations, Eqs. (|2T))) , for the mean 
fields m z A and rn^ x ■ 



Next, we show the case of Nq = 3x3, although we did not use this size of cluster for the scaling analysis in 
this paper. In this case, we have to treat the two clusters given in Fig. I15I1 F and the two corresponding cluster 
Hamiltonians Hci an d Hc 2 ■ The effective fields for the upper cluster in Fig. IT5T II) are given by 



K%%7 = ~2Vim| - 3F 2 m A , h z 2 ' c 4 s 6 8 = -Vim% - 2U 2 m|, h z g ' cS = 0, 
^1,3,5,7 = 4Jm|, h%'f 6!i = 2Jm\, /ig ,cff = 



(C3) 



18 



for the V1-V2 model and 



t z,eff 
"■1,3,5,7 



t z,eff 
"■2,4,6,8 



7 ,cff 



"■1,3,5,7 — tJ "2,4,6 




(C4) 



for the Vdip model. The expressions for the lower cluster 
in Fig. 115( 11) can be obtained by exchanging m^ x and 
m% x in Eqs. JC3J and (iCH) . 

Appendix D: Accuracy estimation of the CMF 
method 

To evaluate the quantitative accuracy of our CMF 
scaling procedure, we compare the results for the sim- 
ple model only with the NN interaction (i.e., the Vy- 
V2 model with V2 = 0) with the corresponding QMC 
data. It is known that no CSS phase appears for V2 = 
and the system exhibits only the direct first-order tran- 
sition from the CS to SF phase i 37 ' 58 Moreover, since 
the quantum fluctuations do not lift the degeneracy at 
(J, h) = (0.5Vi,0) (called the Heisenberg point), the 
value of J c = 0.5Vi docs not change with the cluster 
size -/Ve- 
in Fig. ITuT ah we see that the region of the CS phase 
gradually gets more narrow with increasing the size of the 



cluster used in the CMF calculations, which means that 
the CS phase is first overestimated in the MF theory and 
then gradually improved by taking into account correla- 
tion effects within the cluster. This kind of systematic 
behavior is also seen in other types of cluster extensions 
of the MF theory^ We perform a linear extrapolation 
toward A = 1 using the three samples of Nc = 2x3, 
3x4, and 4x4, which is indicated as the thick solid 
curve in Fig. ITuT ah We can see that the scaled value and 
the QMC data are in good accordance, although the ex- 
trapolation still slightly overestimates the CS phase. The 
linear fits of the three points are quite good as shown in 
Fig. ITuT bh However, it gets worse for very small val- 
ues of J/Vi, which is attributed to the fact that the 
cluster- "shape" dependence becomes more pronounced. 
The scaled values should be improved by using the re- 
sults of larger-size clusters as sample data for the extrap- 
olation. 



* Electronic address: d-yamamoto@rik en.jp| 

1 H. Matsuda and T. Tsuneto, Suppl. Prog. Theor. Phys. 
46, 411 (1970). 

2 K.-S. Liu and M. E. Fisher, J. Low Temp. Phys. 10, 655 
(1972). 

3 T. Matsubara and H. Matsuda, Prog. Theor. Phys. 16, 569 
(1956). 

4 A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. 
Pfau, Phys. Rev. Lett. 94, 160401 (2005). 

5 M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys. 
Rev. Lett. 107, 190401 (2011). 

6 K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. 
Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401 
(2012). 

7 K.-K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe'er, B. 
Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. 
Jin, and J. Ye, Science 322, 231 (2008). 

8 K. Aikawa, D. Akamatsu, M. Hayashi, K. Oasa, J. 
Kobayashi, P. Naidon, T. Kishimoto, M. Ueda, and S. In- 
ouye, Phys. Rev. Lett. 105, 203001 (2010). 

9 J. Deiglmayr, A. Grochola, M. Repp, K. Mortlbauer, C. 



Gluck, J. Lange, O. Dulieu, R. Wester, and M. Wei- 
demiiller, Phys. Rev. Lett. 101, 133004 (2008). 

10 G. G. Batrouni, F. Hebert, and R. T. Scalettar, Phys. Rev. 
Lett. 97, 087209 (2006). 

11 F. J. Burnell, M. M. Parish, N. R. Cooper, and S. L. 
Sondhi, Phys. Rev. B 80, 174519 (2009). 

12 K. Goral, L. Santos, and M. Lewenstein, Phys. Rev. Lett. 
88, 170406 (2002). 

13 D. L. Kovrizhin, G. V. Pai, and S. Sinha, Europhys. Lett. 
72 162 (2005). 

14 P. Sengupta, L. P. Pryadko, F. Alet, M. Troyer, and G. 
Schmid, Phys. Rev. Lett. 94, 207202 (2005). 

15 V. W. Scarola and S. Das Sarma, Phys. Rev. Lett. 95, 
033003 (2005). 

16 S. Yi, T. Li, and C. P. Sun, Phys. Rev. Lett. 98, 260405 
(2007). 

17 I. Danshita and C. A. R. Sade Melo, Phys. Rev. Lett. 103, 
225301 (2009). 

18 B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. 
Zoller, and G. Pupillo, Phys. Rev. Lett. 104, 125301 
(2010). 



19 



(a) 



5 
^5 







26 






27 






28 


: — mf 

CMF-2x2 




29 






30 


A-»l 




31 


o QMC 




32 






33 






34 



0.20 0.25 0.30 0.35 0.40 0.45 0.50 



J IV, 



(b) 





-0.5 




-0.6 » 


o 




\ 


-0.7 ; 


>. o 




-s: 




i 


-0.8 \ 




-0.9 i 




0.0 



FIG. 16: (a) The CS-SF phase boundaries for the NN inter- 
action model by the MF (thin solid line) and CMF methods 
with Nc = 2x2 (dashed line) and 4x4 (dash-dotted line). 
The thick solid line indicates the scaled value of the CMF 
data with the three largest clusters of Nc = 2 x 3, 3 x 4, 
and 4x4. For comparison, we also plot the numerical data 
obtained by using the QMC method based on the directed 
loop algorithm f Ref. [59h for 16 x 16 lattice sites (circles). The 
value J c = 0.5 Vi is constant with respect to Nc as well as 
ho = 2Vi. (b) Cluster-size scalings (A — > 1) of the CMF data 
for the values of h/ho at the CS-SF first-order transition of 
the NN interaction model. We plot the values on the lower 
branch of the transition lines h — ±h' c at J/J c = 0.8, 0.6, and 
0.4. The solid lines are the linear fits of the three samples 
(No = 2 x 3, 3 x 4, and 4 x 4). 
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